pro hisa_spectra, bgpsnum=bgpsnum,hardcopy=hardcopy



;cloudnum=1000
;glon=[32.975,32.975]
;glat=[-0.0714,-0.0714]
;vlsr=[89.0,89.0]

;restore,'/data13/nordhaus/gbt_data/ammonia/analysis/l033/l033_nh3_flag_1_nodup.sav';

restore,'/data13/nordhaus/bolocam_analysis/innergalnh3/results/innergalnh3.sav'

;read in bgpsnum to start at, if specified.
if n_e(bgpsnum) gt 0 then begin
    wbgpsnum = where(innergalnh3.cloudnum ge bgpsnum)
    innergalnh3 = innergalnh3[wbgpsnum]
endif

;no CO file for l009 field, so skip them for this instant
w=where(innergalnh3.glon gt 15.)
innergalnh3_nol009=innergalnh3[w]

b=sort(innergalnh3_nol009.cloudnum)
cloudnum=innergalnh3_nol009[b].cloudnum
glon=innergalnh3_nol009[b].glon ;positions of nh3 pointings
glat=innergalnh3_nol009[b].glat ;positions of nh3 pointings
vlsr=innergalnh3_nol009[b].vlsr


num=n_e(glon)
;;for now define fake cloudnum array
;cloudnum=replicate(1000,num)

moleclinefilename=strarr(num)
hifilename=strarr(num)
vgpspath='/data13/nordhaus/VGPS/'
sgpspath='/data13/nordhaus/SGPS/'
grspath='/data13/nordhaus/GRS/'
for i=0,num-1 do begin
    ;something here for l009 field.
    if glon[i] ge 19.0 and glon[i] le 20.5 then moleclinefilename[i]=grspath+'grs-20-cube.fits'
    if glon[i] gt 20.5 and glon[i] le 21.5 then moleclinefilename[i]=grspath+'grs-21-cube.fits'
    if glon[i] gt 21.5 and glon[i] le 23.0 then moleclinefilename[i]=grspath+'grs-22-cube.fits'
    if glon[i] ge 31.0 and glon[i] le 32.5 then moleclinefilename[i]=grspath+'grs-32-cube.fits'
    if glon[i] gt 32.5 and glon[i] le 33.5 then moleclinefilename[i]=grspath+'grs-33-cube.fits'
    if glon[i] gt 33.5 and glon[i] le 35.0 then moleclinefilename[i]=grspath+'grs-34-cube.fits'
    if glon[i] ge 51.0 and glon[i] le 52.5 then moleclinefilename[i]=grspath+'grs-52-cube.fits'
    if glon[i] gt 52.5 and glon[i] le 53.5 then moleclinefilename[i]=grspath+'grs-53-cube.fits'
    if glon[i] gt 53.5 and glon[i] le 54.5 then moleclinefilename[i]=grspath+'grs-54-cube.fits'
    if glon[i] gt 54.5 and glon[i] le 56.0 then moleclinefilename[i]=grspath+'grs-55-cube.fits'

    if glon[i] ge 7.3 and glon[i] le 10.2 then hifilename[i]=sgpspath+'G8.750.0.0.HI.sgps.fits'
    if glon[i] ge 18.5 and glon[i] le 23.5 then hifilename[i]=vgpspath+'MOS_021.Tb.fits'
    if glon[i] gt 30.5 and glon[i] le 35.5 then hifilename[i]=vgpspath+'MOS_033.Tb.fits'
    if glon[i] gt 50.5 and glon[i] le 55.0 then hifilename[i]=vgpspath+'MOS_053.Tb.fits'
    if glon[i] gt 55.0 and glon[i] le 59.5 then hifilename[i]=vgpspath+'MOS_057.Tb.fits'
end


crap=strarr(1)

oldhifilename='junk'
oldmoleclinefilename='junk'
oldcloudnum=0
cloudnumcounter=0
;loop over all sources
for i=0,num-1 do begin
    
    ;check if we've already read in the fits file
    if hifilename[i] ne oldhifilename then begin
        ;read in HI
        hi=readfits(hifilename[i],hihdr)

        hicrpix1 = sxpar(hihdr,'CRPIX1')-1. ;-1 to correct for difference between numbering in header vs idl
        hicrpix2 = sxpar(hihdr,'CRPIX2')-1.
        hicrpix3 = sxpar(hihdr,'CRPIX3')-1.
        hicrval1 = sxpar(hihdr,'CRVAL1')
        hicrval2 = sxpar(hihdr,'CRVAL2')
        hicrval3 = sxpar(hihdr,'CRVAL3')
        hicdelt1 = sxpar(hihdr,'CDELT1')
        hicdelt2 = sxpar(hihdr,'CDELT2')
        hicdelt3 = sxpar(hihdr,'CDELT3')
        hinaxis1 = sxpar(hihdr,'NAXIS1')
        hinaxis2 = sxpar(hihdr,'NAXIS2')
        hinaxis3 = sxpar(hihdr,'NAXIS3')

        hix = lindgen(hinaxis1)
        hiy = lindgen(hinaxis2)
        hiz = lindgen(hinaxis3)
        hil = (hix-hicrpix1)*hicdelt1+hicrval1
        hib = (hiy-hicrpix2)*hicdelt2+hicrval2
        hiv = (hiz-hicrpix3)*hicdelt3+hicrval3

        extast,hihdr,hiastr
    endif
    ;convert glon glat to pixel coordinates for HI image
    ad2xy,glon[i],glat[i],hiastr,hiw,hiq
    ;pull out HI spectrum
    hivel=hi[hiw,hiq,*,0]

    ;check if we've already read in the fits file
    if moleclinefilename[i] ne oldmoleclinefilename then begin        
        ;read in MOLECLINE
        molecline=readfits(moleclinefilename[i],moleclinehdr)

        moleclinecrpix1 = sxpar(moleclinehdr,'CRPIX1')-1. ;-1 to correct for difference between numbering in header vs idl
        moleclinecrpix2 = sxpar(moleclinehdr,'CRPIX2')-1. ;-1 to correct for difference between numbering in header vs idl
        moleclinecrpix3 = sxpar(moleclinehdr,'CRPIX3')-1. ;-1 to correct for difference between numbering in header vs idl
        moleclinecrval1 = sxpar(moleclinehdr,'CRVAL1')
        moleclinecrval2 = sxpar(moleclinehdr,'CRVAL2')
        moleclinecrval3 = sxpar(moleclinehdr,'CRVAL3')
        moleclinecdelt1 = sxpar(moleclinehdr,'CDELT1')
        moleclinecdelt2 = sxpar(moleclinehdr,'CDELT2')
        moleclinecdelt3 = sxpar(moleclinehdr,'CDELT3')
        moleclinenaxis1 = sxpar(moleclinehdr,'NAXIS1')
        moleclinenaxis2 = sxpar(moleclinehdr,'NAXIS2')
        moleclinenaxis3 = sxpar(moleclinehdr,'NAXIS3')

        moleclinex = lindgen(moleclinenaxis1)
        molecliney = lindgen(moleclinenaxis2)
        moleclinez = lindgen(moleclinenaxis3)
        moleclinel = (moleclinex-moleclinecrpix1)*moleclinecdelt1+moleclinecrval1
        moleclineb = (molecliney-moleclinecrpix2)*moleclinecdelt2+moleclinecrval2
        moleclinev = (moleclinez-moleclinecrpix3)*moleclinecdelt3+moleclinecrval3

        extast,moleclinehdr,moleclineastr
    endif
    ;convert glon glat to pixel coordinates
    ad2xy,glon[i],glat[i],moleclineastr,moleclinew,moleclineq
    ;pull out molecline spectra
    moleclinevel=molecline[moleclinew,moleclineq,*,0]

    ;;scale 13CO to peak of HI emission
    ;scale=max(hivel)/max(moleclinevel)
    ;moleclinevelscaled=moleclinevel*scale

    ;scale HI emission to 13CO peak
    scale2=max(moleclinevel)/max(hivel)
    hivelscaled=hivel*scale2

    ;setup filename in order to save multiple spectra per bgps source
    if cloudnum[i] eq oldcloudnum then begin
        cloudnumcounter = cloudnumcounter+1
        pltfile='spectra/HI_spectrum_'+strtrim(string(cloudnum[i]),2)+'_'+strtrim(string(cloudnumcounter),2)+'.ps'
    endif else begin
        cloudnumcounter = 0
        pltfile='spectra/HI_spectrum_'+strtrim(string(cloudnum[i]),2)+'.ps'
    endelse

    if keyword_set(hardcopy) then begin
        !p.multi=[0,1,0,0,0]
        !x.thick=3.0
        !y.thick=3.0
        !p.thick=3.0
        !p.charsize=1.5
        !p.charthick=3
        set_plot,'ps'
        device,bits_per_pixel=8,yoffset=4.5,xoffset=1.5,/portrait,filename=pltfile,xsize=20,ysize=11,/encapsulated,/color
        ;device,bits_per_pixel=8,yoffset=4.5,xoffset=1.5,/portrait,filename=pltfile,xsize=8,ysize=8,/inches,/encapsulated
    endif

    cgLoadct,39,/silent    

    ;plot this spectra
    plot,moleclinev/1.d3,moleclinevel,xrange=[min(moleclinev/1.d3),max(moleclinev/1.d3)],yrange=[min(moleclinevel)-min(moleclinevel)*0.1,max(moleclinevel)+max(moleclinevel)*0.1],ytitle='Scaled Intensity [K]',xtitle='velocity [km/s]',ystyle=1,xstyle=1,thick=3,line=0,color=75
    oplot,hiv/1.d3,hivelscaled,line=5,thick=3,color=150
    oplot,[vlsr[i],vlsr[i]],[min(moleclinevel)-min(moleclinevel)*0.1,max(moleclinevel)+max(moleclinevel)*0.1],thick=4,line=1
    xyouts,0.78*max(moleclinev/1.d3),0.95*max(moleclinevel),cloudnum[i],charsize=1.2
    plot,[0],[0],/nodata,xrange=[min(moleclinev/1.d3),max(moleclinev/1.d3)],yrange=[min(moleclinevel)-min(moleclinevel)*0.1,max(moleclinevel)+max(moleclinevel)*0.1],ytitle='Scaled Intensity [K]',xtitle='velocity [km/s]',ystyle=1,xstyle=1,thick=3,line=0,/noerase
    if cloudnumcounter ne 0 then xyouts,0.83*max(moleclinev/1.d3),0.89*max(moleclinevel),cloudnumcounter,charsize=1.2

    ;move current filenames to oldfilenames
    oldhifilename=hifilename[i]
    oldmoleclinefilename=moleclinefilename[i]
    oldcloudnum = cloudnum[i]

    ;set up to tab through spectra only if printing to screen
    ;and close file if doing hardcopy
    if keyword_set(hardcopy) eq 0 then begin
        print,'n for next; q for quit'
        read,crap
        if crap eq 'q' then goto,quit
        if crap eq 'n' then goto,next
        next:                   
    endif else begin
        device,/close_file
    endelse


endfor

quit:print,'quitting....'



end
